// Radiosity handling code
// Steps:
//    Sample generation
//        Patches generate radiosity samples for each polygon
//        Trimeshes generate radiosity samples for each vertex
//        Planar surfaces are subdivided into small samples
//    Sample collection
//        Radiosity samples are traced to provide an initial light
//        contribution
//    Transfer creation
//        All radiosity samples are tested against each other.  Transfers
//        are created for any samples capable of diffusing to other samples.
//    Bounce
//        Light is bounced across transfers
//    Finalization
//        Light is stored on to the base lightmaps

#include "cmdlib.h"
#include "mathlib.h"
#include "imagelib.h"
#include "threads.h"
#include "mutex.h"
#include "scriplib.h"

#include "qfiles.h"
#include "shaders.h"
#include "mesh.h"
#include "bspfile.h"

#include "rad.h"
#include "light.h"
#include "lightcommon.h"
#include "../code/qcommon/qfiles.h"

#define	ON_EPSILON	0.1
#define PLANE_EPSILON 0.9999

#define SQR_TWO			1.4142135623730950488016887242097
#define SQR_THREE		1.7320508075688772935274463415059
#define SQR_SIX			2.4494897427831780981972840747059

float basisVectors[] = {
	0.0f,              0.0f,        1.0f,
	-1.0/SQR_SIX,     -1.0/SQR_TWO, 1.0/SQR_THREE,
	-1.0/SQR_SIX,      1.0/SQR_TWO, 1.0/SQR_THREE,
	SQR_TWO/SQR_THREE, 0.0,         1.0/SQR_THREE,
};

static dplane_t patchPlane;
static float radSubdiv = 64.0f;
static float radClusterDist = 64.0f;

static int clusterBytes;

typedef struct
{
	vec3_t dir[MAP_BASIS_VECTORS+1];
} lightSample_t;

typedef struct
{
	vec3_t origin;
	vec3_t normal;

	vec3_t svec, tvec;

	lightSample_t light;

	dplane_t *plane;

	// reportedOrigin is what the sample reports as for triangulation tests
	// For patches, it's an extrapolation of the lightmap coordinates
	// For world surfaces, it's the origin
	vec3_t reportedOrigin;

	vec_t area;

	int cachePlaneIndex;
	int cluster;
} radSample_t;

typedef struct
{
	vec3_t normal;
	float dist;
} cachePlane_t;

static radSample_t *radSamples = NULL;
static int numRadSamples = 0;
static int numTrisurfInputs = 0;
static int numPatchInputs = 0;

static cachePlane_t *cachePlanes = NULL;
static int numCachePlanes = 0;

static radSample_t *AllocRadSample(void)
{
	radSample_t *rs;

	if((numRadSamples % RADSAMPLES_BLOCKSIZE) == 0)
	{
		radSamples = realloc(radSamples, sizeof(radSample_t) * (numRadSamples + RADSAMPLES_BLOCKSIZE));
		if(!radSamples)
			Error("Could not allocate radiosity samples");
	}

	rs = radSamples + numRadSamples;
	numRadSamples++;
	return rs;
}

static cachePlane_t *AllocCachePlane(void)
{
	cachePlane_t *cp;

	if((numCachePlanes % RADSAMPLES_BLOCKSIZE) == 0)
	{
		cachePlanes = realloc(cachePlanes, sizeof(cachePlane_t) * (numCachePlanes + RADSAMPLES_BLOCKSIZE));
		if(!cachePlanes)
			Error("Could not allocate plane cache");
	}

	cp = cachePlanes + numCachePlanes;
	numCachePlanes++;
	return cp;
}

static int CachePlaneIndex(vec3_t normal, float dist)
{
	int i;
	cachePlane_t *cp;

	for(i=0;i<numCachePlanes;i++)
	{
		if(fabs(dist - cachePlanes[i].dist) < ON_EPSILON && DotProduct(normal, cachePlanes[i].normal) > PLANE_EPSILON)
			return i;
	}

	cp = AllocCachePlane();
	VectorCopy(normal, cp->normal);
	cp->dist = dist;

	return cp - cachePlanes;
}


// The way this works: Each point creates a cell enclosing the area where it
// is closer than any other point.  The cell points are interpolated between
// all other cells which share it.  A trifan is then created starting at the
// point itself for interpolation.
#define	MAX_TRI_POINTS		2048

typedef struct
{
	lightSample_t light;
	vec3_t origin;

	int contributedPoints;
} sharedPoint_t;

typedef struct
{
	int			numCells;
	int			numSharedPoints;
	vec3_t		normal;

	radSample_t	*cells[MAX_TRI_POINTS];
	sharedPoint_t sharedPoints;
} triangulation_t;

/*
===============
AllocTriangulation
===============
*/
triangulation_t	*AllocTriangulation (vec3_t normal)
{
	triangulation_t	*t;

	t = malloc(sizeof(triangulation_t));
	t->numCells = 0;
	t->numSharedPoints = 0;

	VectorCopy(normal, t->normal);

//	memset (t->edgematrix, 0, sizeof(t->edgematrix));

	return t;
}

/*
===============
FreeTriangulation
===============
*/
void FreeTriangulation (triangulation_t *tr)
{
	free (tr);
}


/*
============
TriangulatePoints
============
*/
void TriangulatePoints (triangulation_t *trian)
{
}

/*
===============
AddPointToTriangulation
===============
*/
void AddPointToTriangulation (radSample_t *patch, triangulation_t *trian)
{
	int			pnum;

	pnum = trian->numCells;
	if (pnum == MAX_TRI_POINTS)
		Error ("trian->numCells == MAX_TRI_POINTS");
	trian->cells[pnum] = patch;
	trian->numCells++;
}

/*
===============
SampleTriangulation
===============
*/
//#define DEBUG_TRIANGULATION
//#define HARSH_LERPING

void SampleTriangulation (vec3_t point, triangulation_t *trian, vec3_t color)
{
	float best, d;
	int j;
	radSample_t *p1, *p0;
	vec3_t v1;


	// search for nearest point
	best = 99999;
	p1 = NULL;
	for (j=0 ; j<trian->numCells ; j++)
	{
		p0 = trian->cells[j];
		VectorSubtract (point, p0->reportedOrigin, v1);
		d = VectorLength (v1);
		if (d < best)
		{
			best = d;
			p1 = p0;
		}
	}

	if (!p1)
		Error ("SampleTriangulation: no points");

	VectorCopy (p1->light.dir[0], color);

#ifdef DEBUG_TRIANGULATION
	if(best > 8.0f)
	{
		color[0] = 0.0f;
		color[1] = 64.0f;
		color[2] = 64.0f;
	}
#endif
}

/*
=================================================================

  SAMPLE GENERATION

=================================================================
*/

static void SampleFromWinding_r(lFacet_t *facet, winding_t *winding, qboolean patch)
{
	int d;
	vec3_t mins, maxs, dimensions;
	vec3_t splitPlane;
	vec_t dist;
	winding_t *o1, *o2;
	radSample_t *sample;
	float lmoffset[2];
	int leafnum;
	vec3_t origin;

	if(!winding->numpoints)
		Error("winding->numpoints == 0");

	// Check dimensions
	WindingBounds(winding, mins, maxs);
	VectorSubtract(maxs, mins, dimensions);

	splitPlane[0] = splitPlane[1] = splitPlane[2] = 0.0;

	for(d=0;d<3;d++)
	{
		if(maxs[d]-mins[d] > radSubdiv)
		{
			// Too big in this dimension, cut it in half
			splitPlane[d] = 1.0;
			dist = (mins[d]+maxs[d]) * 0.5f;

			o1 = o2 = NULL;
			ClipWindingEpsilon(winding, splitPlane, dist, ON_EPSILON, &o1, &o2);

			// Create samples from both sides
			if(o1)
			{
				if(o1->numpoints == 0)
					Error("frontside numpoints == 0");
				SampleFromWinding_r(facet, o1, patch);
				FreeWinding(o1);
			}
			else
				Error("Split failed?");
			if(o2)
			{
				if(o2->numpoints == 0)
					Error("backside numpoints == 0");
				SampleFromWinding_r(facet, o2, patch);
				FreeWinding(o2);
			}
			else
				Error("Split failed?");

			return;
		}
	}

	// The patch is small enough, convert it into a sample
	WindingCenter(winding, origin);

	leafnum = VL_LightLeafnum(origin);
	// Sample wound up inside a solid, ignore it
	if(dleafs[leafnum].cluster == -1)
		return;

	sample = AllocRadSample();
	memset(sample, 0, sizeof(*sample));

	VectorCopy(origin, sample->origin);
	sample->cluster = dleafs[leafnum].cluster;

	sample->area = WindingArea(winding);
	VectorCopy(facet->plane.normal, sample->normal);

	VectorCopy(facet->textureMatrix[0], sample->svec);
	VectorCopy(facet->textureMatrix[1], sample->tvec);

	// Flatten the S and T vectors against the normal
	dist = DotProduct(sample->svec, sample->normal);
	VectorMA(sample->svec, -dist, sample->normal, sample->svec);
	if(!VectorNormalize(sample->svec, sample->svec))
	{
		printf("WARNING: Degenerate S vector (%4.4f %4.4f %4.4f)\n", sample->origin[0], sample->origin[1], sample->origin[2]);
		printf("%4.4f %4.4f %4.4f\n", facet->textureMatrix[0][0], facet->textureMatrix[0][1], facet->textureMatrix[0][2]);
		printf("%i points on winding\n", winding->numpoints);
	}
	dist = DotProduct(sample->tvec, sample->normal);
	VectorMA(sample->tvec, -dist, sample->normal, sample->tvec);
	if(!VectorNormalize(sample->tvec, sample->tvec))
		printf("WARNING: Degenerate T vector (%4.4f %4.4f %4.4f)\n", sample->origin[0], sample->origin[1], sample->origin[2]);

	// Flatten the T vector against the S vector
	dist = DotProduct(sample->tvec, sample->svec);
	VectorMA(sample->tvec, -dist, sample->svec, sample->tvec);
	if(!VectorNormalize(sample->tvec, sample->tvec))
		printf("WARNING: S/T vectors parallel\n");

	if(patch)
	{
		lmoffset[0] = facet->lightmapCoords[1][0] - DotProduct(facet->points[1], facet->lightmapMatrix[0]);
		lmoffset[1] = facet->lightmapCoords[1][1] - DotProduct(facet->points[1], facet->lightmapMatrix[1]);

		sample->reportedOrigin[0] = (DotProduct(sample->origin, facet->lightmapMatrix[0]) + lmoffset[0]) * RAD_PATCH_REPORT_SCALE;
		sample->reportedOrigin[1] = (DotProduct(sample->origin, facet->lightmapMatrix[1]) + lmoffset[1]) * RAD_PATCH_REPORT_SCALE;
		sample->reportedOrigin[2] = 0.0;
		sample->cachePlaneIndex = -1;
	}
	else
	{
		sample->cachePlaneIndex = facet->planeIndex;
		VectorCopy(sample->origin, sample->reportedOrigin);
	}
}

// SampleFromFacet - Makes samples from triangle on a given facet
static void SampleFromFacetTriangle(lFacet_t *facet, vec3_t v1, vec3_t v2, vec3_t v3, qboolean patch)
{
	winding_t *winding;
	int d;

	// Create a base winding
	winding = AllocWinding(3);
	winding->numpoints = 3;
	VectorCopy(v1, winding->p[0]);
	VectorCopy(v2, winding->p[1]);
	VectorCopy(v3, winding->p[2]);

	// Create samples from it
	SampleFromWinding_r(facet, winding, patch);
	FreeWinding(winding);
}

// RadSamplePatchFacets - Generates facets from 
void RadSamplePatchFacets(lFacet_t *facets, int numFacets)
{
	numPatchInputs += numFacets;
	while(numFacets)
	{
		SampleFromFacetTriangle(facets, facets->points[0], facets->points[1], facets->points[2], qtrue);
		if(facets->numpoints == 4)
			SampleFromFacetTriangle(facets, facets->points[0], facets->points[3], facets->points[2], qtrue);
		facets++;
		numFacets--;
	}

}

int RadPlaneCache(lFacet_t *facet)
{
	int cpi;
	int cptested;

	vec_t actualDist;

	actualDist = DotProduct(facet->points[0], facet->plane.normal);

	cpi = CachePlaneIndex(facet->plane.normal, facet->plane.dist);
	cptested = CachePlaneIndex(facet->plane.normal, actualDist);

	if(cpi != cptested)
		Error("cpi != cptested");

	return cpi;
}

void RadSampleTrisurfFacets(lFacet_t *facets, int numFacets, int planenum)
{
	numTrisurfInputs += numFacets;
	while(numFacets)
	{
		facets->planeIndex = planenum;
		SampleFromFacetTriangle(facets, facets->points[0], facets->points[1], facets->points[2], qfalse);
		if(facets->numpoints == 4)
			SampleFromFacetTriangle(facets, facets->points[0], facets->points[3], facets->points[2], qfalse);
		facets++;
		numFacets--;
	}
}

void RadInitialLight(int sampleNum)
{
	radSample_t *rs;
	light_t *light;
	vec3_t dir;
	float factor;
	traceWork_t tw;
	vec3_t sample;
	vec3_t finalSample;

	tw.patchshadows = qtrue;

	VectorClear(finalSample);

	rs = radSamples + sampleNum;

	for ( light = lights ; light ; light = light->next )
	{
		// Get the direction to the light
		VectorSubtract(light->origin, rs->origin, dir);

		// If the light's behind the sample, ignore it
		if(DotProduct(dir, rs->normal) < 0.0)
			continue;

		VectorNormalize(dir, dir);
		factor = DotProduct(dir, rs->normal);

		// Get the light's contribution
		if (!LightContributionToPoint(light, rs->origin, sample, &tw ))
			continue;

		// Add it to the sample result
		VectorMA(finalSample, factor, sample, finalSample);
	}

	VectorCopy(finalSample, rs->light.dir[0]);
}

static qboolean CullPVS(int fromCluster, int toCluster)
{
	unsigned char *pvs;

	if(!numVisBytes)
		return qtrue;

	pvs = visBytes + 8 + (fromCluster * clusterBytes);
	if( !(pvs[toCluster>>3] & (1<<(toCluster&7))) )
		return qfalse;

	return qtrue;
}

static qboolean CullTransfer(radSample_t *rs, radSample_t *toRS)
{
	vec3_t vec;
	vec3_t midpoint;
	int leafnum;

	VectorSubtract(toRS->origin, rs->origin, vec);

	// Don't transfer if the source is behind the destination sample
	if(DotProduct(vec, toRS->normal) > 0.0)
		return qfalse;

	// Don't transfer if the destination is behind the source sample
	if(DotProduct(vec, rs->normal) < 0.0)
		return qfalse;

	// Check PVS continuity
	if(!CullPVS(rs->cluster, toRS->cluster))
		return qfalse;

	// See if the halfway point is in a solid
	VectorScale(rs->origin, 0.5, midpoint);
	VectorMA(midpoint, 0.5, toRS->origin, midpoint);

	leafnum = VL_LightLeafnum(midpoint);
	if(dleafs[leafnum].cluster == -1)
		return qfalse;

	// See if either of the two quarter points are in a solid
	VectorScale(rs->origin, 0.25, midpoint);
	VectorMA(midpoint, 0.75, toRS->origin, midpoint);

	leafnum = VL_LightLeafnum(midpoint);
	if(dleafs[leafnum].cluster == -1)
		return qfalse;

	VectorScale(rs->origin, 0.75, midpoint);
	VectorMA(midpoint, 0.25, toRS->origin, midpoint);

	leafnum = VL_LightLeafnum(midpoint);
	if(dleafs[leafnum].cluster == -1)
		return qfalse;

	return qtrue;
}

void RadMakeTransfers(int sampleNum)
{
	int i, leafnum;
	trace_t trace;
	traceWork_t tw;
	radSample_t *rs;
	radSample_t *toRS;
	float dist1, dist2;
	vec3_t vec;
	int successful = 0;
	unsigned char *pvs;
	vec3_t midpoint;

	tw.patchshadows = 1;

	rs = radSamples + sampleNum;

	dist1 = DotProduct(rs->normal, rs->origin);

	for(i=0;i<numRadSamples;i++)
	{
		toRS = radSamples + i;

		if(!CullTransfer(rs, toRS))
			continue;

		// Trace it
		//TraceLine(rs->origin, toRS->origin, &trace, qfalse, &tw);

		if(!trace.passSolid)
			successful++;
	}
}


int RadNumSamples(void)
{
	return numRadSamples;
}

void *RadTriangulateSequential(int firstSample, int numSamples, vec3_t normal)
{
	radSample_t *rs;
	triangulation_t *tri;
	int i;

	rs = radSamples + firstSample;

	tri = AllocTriangulation(normal);

	// Add each point
	while(numSamples)
	{
		AddPointToTriangulation(rs, tri);
		rs++;
		numSamples--;
	}

	// Triangulate
	TriangulatePoints(tri);

	if(!tri->numCells)
	{
		FreeTriangulation(tri);
		tri = NULL;
	}

	return tri;
}

static qboolean PointNearSequence(radSample_t *baseSample, int firstSample, int numSamples)
{
	vec3_t v;
	radSample_t *test;

	int i;

	test = radSamples + firstSample;

	while(numSamples)
	{
		VectorSubtract(test->origin, baseSample->origin, v);
		if(VectorLength(v) < radClusterDist)
			return qtrue;

		test++;
		numSamples--;
	}

	return qfalse;
}

void *RadTriangulatePlane(int firstSample, int numSamples, int planeindex)
{
	radSample_t *rs;
	triangulation_t *tri;
	int i;

	tri = AllocTriangulation(cachePlanes[planeindex].normal);

	// Add any point with this planeindex
	for(i=0;i<numRadSamples;i++)
	{
		rs = radSamples + i;
		if(rs->cachePlaneIndex != planeindex)
			continue;

		// Reject anything not near this point cluster to avoid overflows
		if(i < firstSample || i >= (firstSample + numSamples))
			if(!PointNearSequence(rs, firstSample, numSamples))
				continue;

		AddPointToTriangulation(rs, tri);
	}

	// Triangulate
	TriangulatePoints(tri);

	if(!tri->numCells)
	{
		FreeTriangulation(tri);
		tri = NULL;
	}

	return tri;
}

void RadTriangulatePoint(void *surf, vec3_t point, vec3_t out)
{
	SampleTriangulation(point, surf, out);
}

void RadFreeTriangulation(void *ptr)
{
	FreeTriangulation(ptr);
}


// CalcRadiosity - Entry point for all radiosity calculations
void CalcRadiosity(void)
{
#ifndef DO_RADIOSITY
	return;
#endif

	printf("=== Radiosity ===\n");
	printf("%6d patch samples\n", numPatchInputs);
	printf("%6d planar samples\n", numTrisurfInputs);
	printf("%6d subdivisions\n", numRadSamples - (numTrisurfInputs+numPatchInputs));
	printf("%6d radiosity samples\n", numRadSamples);
	printf("%6d cached planes\n", numCachePlanes);

	// Set up the fake patch plane
	patchPlane.dist = 0.0f;
	patchPlane.normal[0] = 0.0f;
	patchPlane.normal[1] = 0.0f;
	patchPlane.normal[2] = 1.0f;

	printf("--- RadMakeTransfers ---\n");
	if(numVisBytes)
		memcpy(&clusterBytes, visBytes + 4, 4);
	else
		printf("WARNING: No vis data, PVS optimization disabled!\n");

	//RunThreadsOnIndividual(numRadSamples, qtrue, RadMakeTransfers);

	printf("--- RadInitialLight ---\n");
	RunThreadsOnIndividual(numRadSamples, qtrue, RadInitialLight);
}

void ApplyRadiosity(void)
{
	_printf("--- ApplyRadiosity ---\n");
	RunThreadsOnIndividual(numDrawSurfaces, qtrue, VL_RadSampleLightmaps);
}
